www.gusucode.com > 2D LINEAR TRUSS FINITE ELEMENT 工具箱matlab源码程序 > 2D LINEAR TRUSS FINITE ELEMENT/Linear_Truss_Finite_Element_2D.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% 2D LINEAR TRUSS FINITE ELEMENT %%% %%% BY: CHRISTOPHER WONG. %%% %%% UNIVERSITY OF NEW HAMPSHIRE %%% %%% crswong888@gmail.com %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% This program computes the displacements, support reactions, and stresses for a 2D truss finite element %%% The formulation is for 2-noded elements with first order linear shape functions %%% The SLE is solved for by eliminating i-th rows and i-the columns for which u(i)=0 (boundary conditions) %%% A user is required to input truss information %%% The algorithm should generally work for any truss structure clear all % clear workspace to avoid conflicts %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% PREPROCESSING (USER INPUT) %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %// Establish coordiantes and global dof's of nodes, node = [name, dof x, x, dof y, y]: node = [ 1, 1, 0, 2, 0 ; 2, 3, -500, 4, 0 ; 3, 5, 400, 6, -300 ] ; % x and y coordinates in mm n_node = size(node,1) ; % establish number of nodes (**NO INPUT**) %// Discretization - assign elements to global nodes, elem = [name, node k, node l]: elem = [ 1, node(2,:), node(1,:) ; 2, node(1,:), node(3,:) ] ; n_elem = size(elem,1) ; % establish number of elements (**NO INPUT**) %// Establish geometric and material properties: E(1:6) = 70 ; % Young's modulus (kN/mm^2) A(1:6) = 200 ; % Cross-sectional area (mm^2) %// Establish applied loads: P = sym('P',[2*n_node 1]) ; % initiate load vector as symbolic array (**NO INPUT**) P(1) = -20 ; % load applied to dof 1 (kN) P(2) = -19 ; % load applied to dof 2 (kN) %// Establish displacement boundary conditions: u = sym('u',[2*n_node 1]) ; % initiate displacement vector as symbolic array (**NO INPUT**) u(3:6,1) = 0 ; % dofs 3 through 6 pinned %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% GENERIC ALGORITHM FOR FEA OF 2D TRUSS %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %// Calculate length and the factors l and m for each element: for i = 1:n_elem L(i) = sqrt( ( elem(i,9) - elem(i,4) )^2 + ( elem(i,11) - elem(i,6) )^2 ) ; % length, mm l(i) = ( elem(i,9) - elem(i,4) ) / L(i) ; % l factor (unitless) m(i) = ( elem(i,11) - elem(i,6) ) / L(i) ; % m factor (unitless) end %// Construct global stiffness matrix from element local matrices: K = zeros(2*n_node,2*n_node); % initiate global stiffness matrix as 2D array of zero's for i = 1:n_elem k = E(i) * A(i) / L(i) * [ l(i)^2 l(i)*m(i) -l(i)^2 -l(i)*m(i) ; l(i)*m(i) m(i)^2 -l(i)*m(i) -m(i)^2 ; -l(i)^2 -l(i)*m(i) l(i)^2 l(i)*m(i) ; -l(i)*m(i) -m(i)^2 l(i)*m(i) m(i)^2 ] ; % element local stiffness matrix (kN/mm) index = [ elem(i,3) elem(i,5) elem(i,8) elem(i,10) ] ; % global indices (dofs) for local stiffness matrix K(index,index) = K(index,index) + k ; % add local stiffness matrices to respective global indices (kN/mm) end %// Establish which i-th rows and i-th columns will be kept for the subsequent calculation of 'u': not_eliminated = false(2*n_node) ; % initiate a false (0) logical array corresponding to the size of K for i = 1:2*n_node if u(i) ~= 0 % condition for which a row or column is NOT eliminated, i.e., kept not_eliminated(i) = true ; % set the i-th logical entry to true (1) end end %// Solve for unknown displacements by using the elimination approach with 'not_eliminated' rows and columns: u(not_eliminated) = inv(K(not_eliminated,not_eliminated)) * P(not_eliminated) ; % displacement solution (mm) u = double(u) ; % convert symbolic array to double precision %// Solve for unknown reaction components of the vector 'P': P = K * u ; % includes applied loads and reactions (kN) P = double(P) ; % convert symbolic array to double precision %// Determine elemental strains and stresses: for i = 1:n_elem epsilon(i) = 1 / L(i) * [ -l(i) -m(i) l(i) m(i) ] * [ u(elem(i,3)) ; u(elem(i,5)) ; u(elem(i,8)) ; u(elem(i,10)) ] ; % element strain (unitless) sigma(i) = E(i) * epsilon(i) ; % element stress (kN/mm^2) end